close; clear; clc;
Nt = 50; Nl = 500; 
nl = 100;
Fai = pi / 2;
u = 4*pi*10^(-7); p1 = 100;
i = sqrt(-1);        
r = 2 * nl;
PE = 1.0 * r;
for it = 1:Nt
    w(it) = 10 ^ (it/10);
    k1(it) = sqrt(-i*w(it)*u/p1);
    
    Ex(it) = PE * p1 / (pi*r^3);
    ExA(it) = abs(Ex(it));
    ExB(it) = angle(Ex(it));
    Hy(it) = -PE / (pi*r^3 * k1(it));
    HyA(it) = abs(Hy(it));
    HyB(it) = angle(Hy(it));
    Hz(it) = 3 * PE * p1 / (2*pi*r^4 * i*w(it)*u);
    HzA(it) = abs(Hz(it));
    HzB(it) = angle(Hz(it));
    pEx(it) = pi*r^3 / PE * Ex(it)/p1;
    pExHy(it) = (ExA(it)/HyA(it))^2/(w(it)*u)/p1;
    pHy(it) = w(it)*u*pi^2*r^6 / PE^2 *HyA(it)^2/p1;
    pHz(it) = w(it)*u*2*pi*r^4 / (3*PE) *HzA(it)/p1;
    
end
figure(1)
loglog(w(:),pHy(:))
   


    